"""
CBMPy: Constraint Based Modelling in Python (http://pysces.sourceforge.net/cbm)
============
Copyright (C) 2009-2012 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>

Author: Brett G. Olivier
Contact email: bgoli@users.sourceforge.net
Last edit: $Author: bgoli $ ($Id: CBXML.py 990 2012-02-15 13:52:11Z bgoli $)

"""

import os, time, numpy
from getpass import getuser
import xml.etree.ElementTree as ELTree
from xml.dom.minidom import getDOMImplementation

if vars().has_key('cDir'):
    cDir=vars()['cDir']
else:
    cDir=os.getcwd()

import CBModel
from CBVersion import __DEBUG__, __version__

_HAVE_SBML_ = None
try:
    import libsbml
    _HAVE_SBML_ = True
except ImportError:
    print 'SBML support not available, install libSMBL with the Python bindings (sbml.org)'
    _HAVE_SBML_ = False

_TEMP_XML_FILE_ = '_tmpxml.tmp'
FBA_NS = 'http://www.sbml.org/sbml/level3/version1/fba/version1'
XHTML_NS = 'http://www.w3.org/1999/xhtml'

global SBML_LEVEL, SBML_VERSION 
SBML_LEVEL = 2
SBML_VERSION = 4

def sbml_readSBML2FBA(fname, work_dir=cDir, return_sbml_model=False):
    """
    Read in an SBML Level 2 file with FBA annotation where and return either a CBM model object
    or a (cbm_mod, sbml_mod) pair if return_sbml_model=True

     - *fname* is the filename
     - *work_dir* is the working directory
     - *return_sbml_model* [default=False] return a a (cbm_mod, sbml_mod) pair

    """

    assert _HAVE_SBML_, "\nSBML not available ... install libSBML with Python bindings for SBML support"
    D = libsbml.readSBML(os.path.join(work_dir,fname))
    M = D.getModel()
    if M == None:
        print 'Invalid SBML no model will be generated'
        if not return_sbml_model:
            return None
        else:
            return None, D        
        
    model_id = M.getId()
    model_name = M.getName()
    model_description = libsbml.XMLNode_convertXMLNodeToString(M.getNotes())

    __HAVE_FBA_ANOT__ = False
    ANOT = M.getAnnotationString()
    if len(ANOT) > 0:
        __HAVE_FBA_ANOT_OBJ__ = True
        __HAVE_FBA_ANOT_BNDS__ = True
        __HAVE_FBA_ANOT__ = True
        if '<fba:listOfConstraints/>' in ANOT:
            __HAVE_FBA_ANOT_BNDS__ = False
        if '<fba:listOfObjectives/>' in ANOT:
            __HAVE_FBA_ANOT_OBJ__ = False    
        if not __HAVE_FBA_ANOT_BNDS__ and not __HAVE_FBA_ANOT_OBJ__:
            __HAVE_FBA_ANOT__ = False        
        F = file(os.path.join(work_dir, _TEMP_XML_FILE_),'w')
        F.write(ANOT)
        F.flush()
        F.close()
    del ANOT

    SPEC = []
    for s in range(M.getNumSpecies()):
        SBSp = M.getSpecies(s)
        id = SBSp.getId()
        boundCon = False
        if SBSp.getBoundaryCondition():
            if __DEBUG__: print 'Real boundary metabolite: %s' % id
            boundCon = True
            # DODGY AS A LAST RESORT
        elif id[-2:] == '_b':
            if __DEBUG__: print 'Fake boundary (_b) metabolite: %s' % id
            boundCon = True
        CF = '' # chemical formula
        NM = SBSp.getName() # get name
        # to strip a BiGG file see CBTools
        S = CBModel.Species(SBSp.getId(), boundary=boundCon, name=NM, value=SBSp.getInitialAmount(), compartment=SBSp.getCompartment(), charge=SBSp.getCharge(), chemFormula=CF)
        S.setAnnotation('note', libsbml.XMLNode_convertXMLNodeToString(SBSp.getNotes()))
        SPEC.append(S)
    boundary_species = [s.getPid() for s in SPEC if s.is_boundary]

        ##  CBModel.Species{id : {'name' : SBSp.getName(),
                          ##  'id' : id,
                          ##  'constant' : SBSp.getConstant(),
                          ##  'boundary' : SBSp.getBoundaryCondition(),
                          ##  'concentration' : SBSp.getInitialConcentration(),
                          ##  'amount' : SBSp.getInitialAmount(),
                          ##  'compartment' : SBSp.getCompartment()
                            ##  }
                        ##  })

    spec_id = [s.getPid() for s in SPEC]

    REAC = []
    reactionIDs = []
    reactionsReversability = []
    for r in range(M.getNumReactions()):
        SBRe = M.getReaction(r)
        R_id = SBRe.getId()
        reagents = []
        EXREAC = False
        reactionIDs.append(R_id)
        for re in range(SBRe.getNumReactants()):
            spec = SBRe.getReactant(re).getSpecies()
            stoi = -SBRe.getReactant(re).getStoichiometry()
            reagents.append((stoi,spec))
            if spec in boundary_species:
                EXREAC = True
        for pr in range(SBRe.getNumProducts()):
            spec2 = SBRe.getProduct(pr).getSpecies()
            stoi2 = SBRe.getProduct(pr).getStoichiometry()
            reagents.append((stoi2,spec2))
            if spec2 in boundary_species:
                EXREAC = True
        R = CBModel.Reaction(SBRe.getId(), SBRe.getName(), reversible=SBRe.getReversible())
        reactionsReversability.append(SBRe.getReversible())
        for r in reagents:
            rgtmp = CBModel.Reagent(SBRe.getId()+r[1], r[1], r[0])
            R.addReagent(rgtmp)
            if R.getPid() not in SPEC[spec_id.index(r[1])].reagent_of:
                SPEC[spec_id.index(r[1])].reagent_of.append(R.getPid())
        if EXREAC:
            R.is_exchange = True
        R.setAnnotation('note', libsbml.XMLNode_convertXMLNodeToString(SBRe.getNotes()))
        REAC.append(R)
        
    CONSTR = []
    OBJFUNC = []
    objfunc = {}
    if __HAVE_FBA_ANOT__:
        root = ELTree.ElementTree(file=os.path.join(work_dir, _TEMP_XML_FILE_))
        root_i = root.getiterator()
        for ri in root_i:
            if ri.tag == '{http://www.sbml.org/sbml/level3/version1/fba/version1}fluxBalance':
                if __DEBUG__: print ri.tag
                rootfba = ELTree.ElementTree(ri)
                root_fba_i = rootfba.getiterator()
        constraints = []
        for ret in root_fba_i:
            if __HAVE_FBA_ANOT_BNDS__ and ret.tag == '{http://www.sbml.org/sbml/level3/version1/fba/version1}listOfConstraints':
                if __DEBUG__: print ret.tag
                chld = ret.getchildren()
                for c in chld:
                    if __DEBUG__: print '\t%s' % c.tag
                    attrib = c.attrib
                    for a in attrib:
                        attrib.update({a.replace('{http://www.sbml.org/sbml/level3/version1/fba/version1}','') :  attrib.pop(a)})
                    constraints.append(attrib)
            # black magic happening here needs some generalisation
            elif __HAVE_FBA_ANOT_OBJ__ and ret.tag == '{http://www.sbml.org/sbml/level3/version1/fba/version1}listOfObjectives':
                if __DEBUG__: print ret.tag
                if len(ret) == 0:
                    print 'No objectives in listOfObjectives'
                else:
                    ##  print 'Adding objectives', ret
                    obj = ret[0]
                    ##  print obj.attrib
                    if obj.attrib.has_key('{http://www.sbml.org/sbml/level3/version1/fba/version1}type'):
                        type = obj.attrib['{http://www.sbml.org/sbml/level3/version1/fba/version1}type']
                    else:
                        type = None
                    ##  raw_input(type)
                    id = obj.attrib['id']
                    objfunc = obj.getchildren()[0].getchildren()[0].attrib
                    ##  raw_input(objfunc)
                    if __DEBUG__: print objfunc
                    for a in objfunc:
                        objfunc.update({a.replace('{http://www.sbml.org/sbml/level3/version1/fba/version1}','') :  objfunc.pop(a)})
                    objfunc.update({'type' : type})
                    objfunc.update({'id' : id})
                    if __DEBUG__: print objfunc
                    objective_function = objfunc
        
        print objfunc
        if len(objfunc.keys()) > 0:
            # THIS IS A HACK FOR A SINGLE OBJECTIVE FUNCTION ... GENERALISE
            OBJFUNC = [CBModel.Objective(objfunc['id'], objfunc['type'])]
            OBJFUNC[0].addFluxObjective(CBModel.FluxObjective(objfunc['reaction']+'_obj',\
                                        objfunc['reaction'], float(objfunc['coefficient'])))
        else:
            OBJFUNC = []

        cntr = 0
        boundReactionIDs = []
        for c in constraints:
            if c.has_key('id'):
                metaid = c['id']
            else:
                O = 'C_%i' % cntr
                if c['operation'] in ['greater','greaterEqual','>','>=']:
                    O = 'lower'
                elif c['operation'] in ['less','lessEqual','<','<=']:
                    O = 'upper'
                elif c['operation'] in ['equal','=']:
                    O = 'equal'
                metaid = '%s_%s_bnd'% (c['reaction'], O)        
            CONSTR.append(CBModel.FluxBound(metaid, c['reaction'],\
                                            c['operation'], float(c['value'])))
            cntr+=1
            if c['reaction'] not in boundReactionIDs:
                boundReactionIDs.append(c['reaction'])

        # undefined flux bounds are given infinite value
        ubcntr = 0
        ##  print boundReactionIDs
        ##  print reactionIDs

        for J in range(len(reactionIDs)):
            ##  print reactionIDs[J], reactionsReversability[J]
            if reactionIDs[J] not in boundReactionIDs:
                metaid = 'UC_%i' % ubcntr
                if reactionsReversability[J]:
                    ##  print 'Adding reversible'
                    CONSTR.append(CBModel.FluxBound(metaid, reactionIDs[J], 'greaterEqual', -numpy.inf))
                else:
                    ##  print 'Adding irreversible'
                    CONSTR.append(CBModel.FluxBound(metaid, reactionIDs[J], 'greaterEqual', 0.0))
                ubcntr += 1 
                metaid = 'UC_%i' % ubcntr                
                CONSTR.append(CBModel.FluxBound(metaid, reactionIDs[J], 'lessEqual', numpy.inf))
                ubcntr += 1

    # build model
    fm = CBModel.Model(model_id)
    fm.name = model_name
    fm.description = model_description
    for s in SPEC:
        fm.addSpecies(s)
    for r in REAC:
        fm.addReaction(r)
    for c in CONSTR:
        fm.addFluxBound(c)
    for o in OBJFUNC:
        fm.addObjective(o, active=True)

    if not __HAVE_FBA_ANOT__:
        print '\nWARNING: No FBA annotations found (no flux bounds or objective functions)!\n'
        time.sleep(1)
    if not return_sbml_model:
        return fm
    else:
        return fm, M

def readSBML2FBA(fname, work_dir=cDir, return_sbml_model=False):
    """
    Read in an SBML Level 2 file with FBA annotation where:

     - *fname* is the filename
     - *work_dir* is the working directory
     - *return_sbml_model* [default=False] return a a (cbm_mod, sbml_mod) pair

    """    
    return sbml_readSBML2FBA(fname, work_dir, return_sbml_model)


def xml_createSBML2FBADoc():
    """
    Create a 'document' to store the SBML2FBA annotation, returns:

     - *DOC* a minidom document

    """
    DOM = getDOMImplementation()
    DOC = DOM.createDocument(FBA_NS, 'fba:fluxBalance', None)
    DOC.documentElement.setAttribute('xmlns:fba', FBA_NS)
    TOP = DOC.documentElement
    LoC = DOC.createElementNS(FBA_NS, 'fba:listOfConstraints')
    TOP.appendChild(LoC)
    LoO = DOC.createElementNS(FBA_NS, 'fba:listOfObjectives')
    TOP.appendChild(LoO)
    return DOC

def xml_viewSBML2FBAXML(document, fname=None):
    """
    Print a minidom XML document to screen or file, arguments:

     - *document* a minidom XML document
     - *fname* [default=None] by default print to screen or write to file fname

    """
    ##  print DOC.toxml('UTF-8')
    ##  print document.toprettyxml(indent=' ',newl='\n',encoding='UTF-8')
    if fname != None:
        FO = file(fname, 'w')
        document.writexml(FO, indent='', addindent='    ', newl='\n', encoding='UTF-8')
        FO.flush()
        FO.close()

def xml_addSBML2FBAFluxBound(document, rid, operator, value, id=None):
    """
    Adds an SBML3FBA flux bound to the document:

     - *document* a minidom XML document created by xml_createSBML2FBADoc
     - *rid* the reaction id
     - *operator* one of ['greater','greaterEqual','less','lessEqual','equal','>','>=','<','<=','=']
     - *value* a float which will be cast to a string using str(value)

    """
    OPER = ['greater','greaterEqual','less','lessEqual','equal','>','>=','<','<=','=']
    assert operator in OPER, '\nBad operator %s should be one of %s' % (operator, str(OPER))
    LoC = document.getElementsByTagName('fba:listOfConstraints')[0]
    F = document.createElementNS(FBA_NS, 'fba:constraint')
    F.setAttributeNS(FBA_NS,'fba:reaction', rid)
    F.setAttributeNS(FBA_NS,'fba:operation', operator)
    value = float(value)
    if numpy.isinf(value):
        if value > 0:
            value = 'inf'
        else:
            value = '-inf'
    else:
        value = str(value)
    F.setAttributeNS(FBA_NS,'fba:value', value)
    if id == None:
        O = 'type'
        if operator in ['greater','greaterEqual','>','>=']:
            O = 'lower'
        elif operator in ['less','lessEqual','<','<=']:
            O = 'upper'
        elif operator in ['equal','=']:
            O = 'equal'
        id = '%s_%s_bnd'% (rid, operator)
    
    F.setAttributeNS(FBA_NS,'metaid', id)
    F.setAttributeNS(FBA_NS,'id', id)
    LoC.appendChild(F)

def xml_createListOfFluxObjectives(document, fluxObjectives):
    """
    Create a list of fluxObjectives to add to an Objective:

     - *document* a minidom XML document created by xml_createSBML2FBADoc
     - *fluxobjs* a list of (rid, coefficient) tuples

    """
    LoF = document.createElementNS(FBA_NS, 'fba:listOfFluxes')
    for F in fluxObjectives:
        FO = document.createElementNS(FBA_NS, 'fba:fluxObjective')
        FO.setAttributeNS(FBA_NS,'fba:reaction', F[0])
        FO.setAttributeNS(FBA_NS,'fba:coefficient', F[1])
        LoF.appendChild(FO)
    return LoF

def xml_createSBL2FBAObjective(document, id, type, fluxObjectives):
    """
    Create a list of fluxObjectives to add to an Objective:

     - *document* a minidom XML document created by xml_createSBML2FBADoc
     - *id* the objective id
     - *type* a string containing the objective sense either: **maximize** or **minimize**
     - *fluxObjectives* a list of (rid, coefficient) tuples

    """
    type = type.lower()
    if type == 'max': type = 'maximize'
    if type == 'min': type = 'minimize'        
    if type in ['maximise', 'minimise']:
        type = type.replace('se','ze')
    assert type in ['maximize', 'minimize'], "\nType must be ['maximize', 'minimize'] not %s" % type
    OBJ = document.createElementNS(FBA_NS, 'fba:objective')
    OBJ.setAttributeNS(FBA_NS,'id', id)
    ##  OBJ.setAttributeNS(FBA_NS,'metaid', id)
    OBJ.setAttributeNS(FBA_NS,'fba:type', type)
    LoF = document.createElementNS(FBA_NS, 'fba:listOfFluxes')
    for F in fluxObjectives:
        FO = document.createElementNS(FBA_NS, 'fba:fluxObjective')
        FO.setAttributeNS(FBA_NS,'fba:reaction', F[0])
        FO.setAttributeNS(FBA_NS,'fba:coefficient', str(F[1]))
        LoF.appendChild(FO)
    OBJ.appendChild(LoF)
    return OBJ

def xml_addSBML2FBAObjective(document, objective, active=True):
    """
    Adds an objective element to the documents listOfObjectives and sets the active attribute:

     - *document* a minidom XML document created by `xml_createSBML2FBADoc`
     - *objective* a minidom XML objective element created with `xml_createSBL2FBAObjective`
     - *active* [default=True] a boolean flag specifiying whether this objective is active

    """

    LoO = document.getElementsByTagName('fba:listOfObjectives')[0]
    if active:
        id = objective.getAttributeNS(FBA_NS,'id')
        LoO.setAttributeNS(FBA_NS,'fba:activeObjective', id)
        print 'Setting active objective: %s' % id
    LoO.appendChild(objective)


def xml_getSBML2FBAannotation(fba, fname=None):
    """
    Takes an FBA model object and returns the SBML3FBA annotation as an XML string:

     - *fba* an fba model object
     - *fname* [default=None] if supplied the XML will be written to file *fname*

    """

    DOC = xml_createSBML2FBADoc()
    for f in fba.flux_bounds:
        ##  print f.value
        xml_addSBML2FBAFluxBound(DOC, f.reaction, f.operation, f.value, f.getPid())

    for o in fba.objectives:
        fluxobjs = [(fo.reaction, fo.coefficient) for fo in  o.fluxObjectives]
        OBJ = xml_createSBL2FBAObjective(DOC, o.id, o.operation, fluxobjs)
        if o.id == fba.objectives[fba.activeObjIdx].id:
            xml_addSBML2FBAObjective(DOC, OBJ, active=True)
        else:
            xml_addSBML2FBAObjective(DOC, OBJ, active=False)
    if fname != None:
        xml_viewSBML2FBAXML(DOC, fname)
        
    return DOC.toprettyxml(indent=' ',newl='\n',encoding='UTF-8')


def sbml_createModel(fba, level=2, version=4):
    """
    Create an SBML model and document:
    
     - *fba* a PySCeSCBM model instance
     - *level* [default=2] SBML level
     - *version* [default=4] SBML version
    
    and returns:
    
     - *model* an SBML model
    
    """
    global SBML_LEVEL,SBML_VERSION
    if level == None or version == None:
        pass
    else:        
        SBML_LEVEL = level
        SBML_VERSION = version
        

    if fba.getPid() == '' or fba.getPid() == None:
        mid0 = 'FBAModel'
    else:
        mid0 = fba.getPid()

    mid = ''
    for l in mid0:
        if l.isalnum():
            mid += l
        else:
            mid += '_'

    if not mid[0].isalpha():
        mid = 'id_' + mid

    mname = fba.getName()
    if  mname in [None, '', ' ']:
        mname = 'SBML_CB_MODEL'

    if libsbml.getLibSBMLVersion() < 40000:
        model = libsbml.Model()
        model.setId(mid)
        model.setName(mname)
        ##  document = libsbml.SBMLDocument()
        ##  document.setLevelAndVersion(SBML_LEVEL, SBML_VERSION)
    else:
        model = libsbml.Model(SBML_LEVEL, version)
        model.setId(mid)
        model.setName(mname)
        
        ##  document = libsbml.SBMLDocument(SBML_LEVEL, SBML_VERSION)
    return model

def sbml_setCompartments(model, compartments):
    if len(compartments.keys()) < 1:
        print 'Warning: no compartments defined adding one called \"cell\"'
        compartments.update({'cell' : {'id' : 'cell', 'dimensions' : 3, 'size' : 1.0}})
    for cs in compartments:
        comp_def = model.createCompartment()
        comp_def.setId(compartments[cs]['id'])
        comp_def.setName(compartments[cs]['id'])
        comp_def.setVolume(float(compartments[cs]['size']))

def sbml_setDescription(model, fba):
    '''
    Sets the model description as a <note> containing `txt` in an HTML paragraph on the model object.
    
     - *model* a libSBML model instance
     - *fba* a PySCeSCBM model instance
    
    '''
    
    ##  try: UseR = getuser()
    ##  except: UseR = ''
    notes = ''
    notes += '<span xmlns="%s">\n' % XHTML_NS
    notes += '<p><span style="font-family: Courier New,Courier,monospace;">%s</span></p>' % fba.description
    notes += '<p><br/><span size="small">Model \"<strong>%s</strong>\" (%s) generated with <a href="http://pysces.sourceforge.net">PySCeS CBM</a> (%s) on %s.</span></p>' % (fba.getPid(),fba.getName(),__version__,time.strftime("%a, %d %b %Y %H:%M:%S"))
    notes += '</span>'
    model.setNotes(notes)

def sbml_setUnits(model, units=None, give_default=False):
    """
    Adds units to the model:
    
    - *model* a libSBML model instance
    - *units* [default=None] a dictionary of units, if None default units are used
    - *give_default* [default=False] if true method returns the default unit dictionary
    
    """
    
    ud = {'area': {0 : {'exponent': 2, 'kind': 'metre', 'multiplier': 1.0, 'scale': 0}},
                'length': {0 : {'exponent': 1, 'kind': 'metre', 'multiplier': 1.0, 'scale': 0}},
                'substance': {0 : {'exponent': 1, 'kind': 'mole', 'multiplier': 1.0, 'scale': 0}},
                'time': {0 : {'exponent': 1, 'kind': 'second', 'multiplier': 1.0, 'scale': 0}},
                'volume': {0 : {'exponent': 1, 'kind': 'litre', 'multiplier': 1.0, 'scale': 0}},
                'mmol_per_gDW_per_hr': {0 : {'exponent': 1, 'kind': 'mole', 'multiplier': 1.0, 'scale': -3},
                                        1 : {'exponent': -1, 'kind': 'gram', 'multiplier': 1.0, 'scale': 0},
                                        2 : {'exponent': -1, 'kind': 'second', 'multiplier': 0.00027777, 'scale': 0}
                                      }
                }
                
    if give_default:
        return ud

    for un in ud.keys():
        vdef = model.createUnitDefinition()
        vdef.setId(un)
        vdef.setName(un)
        for u in range(len(ud[un].keys())):
            vu = model.createUnit()
            vu.setKind(libsbml.UnitKind_forName(ud[un][u]['kind']))
            vu.setMultiplier(ud[un][u]['multiplier'])
            vu.setScale(int(ud[un][u]['scale']))
            vu.setExponent(int(ud[un][u]['exponent']))
            vu.setOffset(0)

def sbml_writeAnnotationsAsBiGGNote_XML(annotations):
    """
    Writes the annotations dictionary as a BiGG compatible SBML <note> with the correct namespace
    so it is XML and SBML L2V4 compatible
    """
    annoSTR = '<body xmlns:html="%s">\n' % XHTML_NS
    for K in annotations.keys():
        val = str(annotations[K]).replace('<','&lt;').replace('>','&gt;')
        ##  val = str(annotations[K]).replace("'",'').replace('<','&lt;').replace('>','&gt;')
        annoSTR += "<html:p>%s: %s</html:p>\n" % (K, val)
    annoSTR += '</body>\n'
    ##  print annoSTR
    return annoSTR


def sbml_setSpecies(model, fba, return_dicts=False):
    """
    Add the species definitions to the SBML object:
    
     - *model* [default=''] a libSBML model instance or can be None if *return_dicts* == True
     - *fba* a PySCeSCBM model instance
     - *return_dicts* [default=False] only returns the compartment and species dictionaries without updated the SBML
     
    returns:
    
     - *compartments* a dictionary of compartments (except when give *return_dicts* argument)
     
    """

    compartments = {}
    species = {}
    for s in fba.species:
        if s.compartment not in compartments.keys():
            if s.compartment == None or s.compartment == '':
                pass
            else:
                compartments.update({s.compartment : {'id' : s.compartment, 'size' : 1.0, 'dimensions' : 3}})
                
        species.update({s.getPid() : {'id' : s.getPid(),
                                       'compartment' : s.compartment,
                                       'name' : s.getName(),
                                       'charge' : s.charge,
                                       'value' : s.value,
                                       'annotation' : s.getAnnotations(),
                                       'boundary' : s.is_boundary,
                                       'chemFormula' : s.chemFormula
                                      }
                        })

    if len(compartments.keys()) < 1:
        print 'Warning: no compartments defined adding one called \"cell\"'
        compartments.update({'cell' : {'id' : 'cell', 'dimensions' : 3, 'size' : 1.0}})

    if return_dicts:
        return compartments, species

    keys=species.keys()
    keys.sort()
    for spe in keys:
        s = model.createSpecies()
        s.setId(species[spe]['id'])
        s.setMetaId(species[spe]['id'])
        s.setName(species[spe]['name'])
        if not (species[spe]['charge'] != None or species[spe]['charge'] != ''):
            s.setCharge(int(species[spe]['charge']))
        if species[spe]['compartment'] == None or species[spe]['compartment'] == '':
            s.setCompartment('cell')
        else:
            s.setCompartment(species[spe]['compartment'])
        if species[spe]['boundary']:
            s.setBoundaryCondition(True)
            #s.setConstant(True)
        else:
            s.setBoundaryCondition(False)
        if species[spe]['value'] == None or species[spe]['value'] == '':
            s.setInitialConcentration(0.0)
            ##  s.setInitialAmount(0.0)
        else:
            s.setInitialConcentration(float(species[spe]['value']))
            ##  s.setInitialAmount(float(species[spe]['value']))
        s.setHasOnlySubstanceUnits(False)
        Ranno = species[spe]['annotation']
        annoSTR = None
        if species[spe]['chemFormula'] != None and species[spe]['chemFormula'] != '':
            annoSTR = sbml_writeAnnotationsAsBiGGNote_XML({'chemFormula' : species[spe]['chemFormula']})
            ##  if SBML_LEVEL >= 2 and SBML_VERSION >= 1:
                ##  annoSTR = sbml_writeAnnotationsAsBiGGNote_XML({'chemFormula' : species[spe]['chemFormula']})
            ##  else:
                ##  pass
        elif Ranno.has_key('chemFormula'):
            annoSTR = sbml_writeAnnotationsAsBiGGNote_XML(Ranno)
            ##  if SBML_LEVEL >= 2 and SBML_VERSION >= 1:
                ##  annoSTR = sbml_writeAnnotationsAsBiGGNote_XML(Ranno)
            ##  else:
                ##  annoSTR = '<span xmlns="%s">\n' % XHTML_NS
                ##  if not (species[spe]['chemFormula'] == None or species[spe]['chemFormula'] == ''):
                    ##  annoSTR += '<p><span>chemFormula</span><span>%s</span></p>' % species[spe]['chemFormula']
                ##  elif Ranno.has_key('chemFormula') and (Ranno['chemFormula'] != None or Ranno['chemFormula'] != ''):
                    ##  annoSTR += '<p><span>chemFormula</span><span>%s</span></p>' % Ranno['chemFormula']
                ##  if annoSTR != '<span xmlns="%s">\n' % XHTML_NS:
                    ##  annoSTR += '</span>'
        if annoSTR != None:
            nres = s.setNotes(annoSTR)
    return compartments


def sbml_setReactions(model, fba, return_dict=False):
    """
    Add the FBA instance reactions to the SBML model

     - *model* an SBML model instance
     - *fba* a PySCeSCBM model instance
     - *return_dict* [default=False] if True do not add reactions to SBML document instead return a dictionary description of the reactions

    """
    
    global SBML_LEVEL, SBML_VERSION
    reactions = {}    

    for r in fba.reactions:
        reactants = []
        products = []
        for rg in r.getStoichiometry():
            if rg[0] > 0.0:
                products.append(rg)
            else:
                reactants.append(rg)
        reactions.update({r.getPid() : {'id' : r.getPid(),
                                         'reactants' : reactants,
                                         'products' : products,
                                         'name' : r.getName(),
                                         'reversible': r.reversible,
                                         'exchange' : r.is_exchange,
                                         'annotation' : r.getAnnotations(),
                                         'compartment' : r.compartment
                                        }
                          })    
    if return_dict:
        return reactions
        
    for rxn in reactions:
        # print 'Adding reaction:', reactions[rxn]['id']
        r = model.createReaction()
        r.setId(reactions[rxn]['id'])
        r.setMetaId(reactions[rxn]['id'])
        r.setName(reactions[rxn]['name'])
        for s in range(len(reactions[rxn]['reactants'])):
            #print '\t' + reactions[rxn]['id'] +' has substrate: ' + s.name + ' (%s)' % abs(rxn.stoichiometry[s.name])
            if libsbml.getLibSBMLVersion() < 40000:
                sref = libsbml.SpeciesReference(reactions[rxn]['reactants'][s][1], abs(float(reactions[rxn]['reactants'][s][0])))
            else:
                sref = libsbml.SpeciesReference(SBML_LEVEL, SBML_VERSION)
                sref.setStoichiometry(abs(float(reactions[rxn]['reactants'][s][0])))
                sref.setSpecies(reactions[rxn]['reactants'][s][1])
            sref.setMetaId('%s_%s' % (reactions[rxn]['id'], reactions[rxn]['reactants'][s][1]))
            r.addReactant(sref)
        for p in range(len(reactions[rxn]['products'])):
            #print '\t' + reactions[rxn]['id'] +' has product: ' + p.name + ' (%s)' % abs(rxn.stoichiometry[p.name])
            if libsbml.getLibSBMLVersion() < 40000:
                pref = libsbml.SpeciesReference(reactions[rxn]['products'][p][1], abs(float(reactions[rxn]['products'][p][0])))
            else:
                pref = libsbml.SpeciesReference(SBML_LEVEL, SBML_VERSION)
                pref.setStoichiometry(abs(float(reactions[rxn]['products'][p][0])))
                pref.setSpecies(reactions[rxn]['products'][p][1])
            pref.setMetaId('%s_%s' % (reactions[rxn]['id'], reactions[rxn]['products'][p][1]))
            r.addProduct(pref)
        Ranno = reactions[rxn]['annotation']
        
        if Ranno.has_key('GENE_ASSOCIATION'):
            Ranno.update({'GENE ASSOCIATION' : Ranno.pop('GENE_ASSOCIATION')})      
        annoSTR = sbml_writeAnnotationsAsBiGGNote_XML(Ranno)            
        ##  if SBML_LEVEL >= 2 and SBML_VERSION >= 1:
            ##  annoSTR = sbml_writeAnnotationsAsBiGGNote_XML(Ranno)
        ##  else:
            ##  annoSTR = '<html xmlns="%s"><head><title/></head><body>\n' % (XHTML_NS)
            ##  annoSTR += '<p>Abbreviation: %s</p>\n' % (reactions[rxn]['id'])
            ##  if Ranno.has_key('Equation'):
                ##  annoSTR += '<p>Equation: %s</p>\n' % (Ranno['Equation'].replace('<', '&lt;').replace('>', '&gt;')) 
            ##  if Ranno.has_key('SUBSYSTEM'):
                ##  annoSTR += '<p>SUBSYSTEM: %s</p>\n' % (Ranno['SUBSYSTEM'])            
            ##  if Ranno.has_key('GENE ASSOCIATION'):
                ##  annoSTR += '<p>GENE ASSOCIATION: %s</p>\n' % (Ranno['GENE ASSOCIATION'])
            ##  elif Ranno.has_key('GENE_ASSOCIATION'):
                ##  annoSTR += '<p>GENE ASSOCIATION: %s</p>\n' % (Ranno['GENE_ASSOCIATION'])            
            ##  annoSTR += '</body></html>\n'
        nres = r.setNotes(annoSTR)

        if reactions[rxn]['reversible']:
            r.setReversible(True)
        else:
            r.setReversible(False)


def sbml_writeSBMLModel(model, filename, directory=None, return_doc=False):
    """
    Writes an SBML model object to file. Note this is an internal SBML method use `sbml_writeSBML2FBA()` to write an FBA model:
    
     - *model* a libSBML model instance
     - *filename* the output filename
     - *directory* [default=None] by default use filename otherwise join, <dir><filename>
     - *return_doc* [default=False] return the SBML document used to write the XML
    
    """
    if directory != None:
        assert os.path.exists(directory), '\n%s does not exist.' % directory
        filename = os.path.join(directory, filename)
    print 'Writing file: %s' % filename

    try: UseR = getuser()
    except: UseR = ''

    global SBML_LEVEL, SBML_VERSION

    if libsbml.getLibSBMLVersion() < 40000:
        document = libsbml.SBMLDocument()
        document.setLevelAndVersion(SBML_LEVEL, SBML_VERSION)
    else:
         document = libsbml.SBMLDocument(SBML_LEVEL, SBML_VERSION)
    ##  if SBML_LEVEL == 2 and SBML_VERSION == 1:
        ##  ns = document.getNamespaces()
        ##  ns.add("http://www.w3.org/1999/xhtml", "html")

    document.setModel(model)
    h1 = '<?xml version="1.0" encoding="utf-8"?>\n'
    h1 += '<!-- Created with PySCeS CBM ('+ __version__ + ') on ' + time.strftime("%a, %d %b %Y %H:%M:%S") + ' by '+UseR+' -->\n'
    F = file(filename, 'w')
    F.write(h1 + document.toSBML())
    F.flush()
    F.close()
    print 'Model exported as: %s' % (filename)

def sbml_writeSBML2FBA(fba, fname, directory=None, sbml_level_version=None):
    """
    Takes an FBA model object and writes it to file as SBML L3 FBA:

     - *fba* an fba model object
     - *fname* the model will be written as XML to *fname*
     - *directory* [default=None] if defined it is prepended to fname
     - *sbml_level_version* [default=None] a tuple containing the SBML level and version e.g. (2,4)

    """
    ##  __DEBUG__ = True
    
    assert _HAVE_SBML_, "\nSBML not available ... install libSBML with Python bindings for SBML support"

    
    if sbml_level_version != None:
        global SBML_LEVEL, SBML_VERSION
        SBML_LEVEL = sbml_level_version[0]
        SBML_VERSION = sbml_level_version[1]
        
    #  create L3 as annotation 
    L3FBA = xml_getSBML2FBAannotation(fba, 'test_xml.xml') # should be None
    
    # create a model    
    SMOD = sbml_createModel(fba, level=SBML_LEVEL, version=SBML_VERSION)
    sbml_setDescription(SMOD, fba)
    sbml_setUnits(SMOD, units=None)
    C = sbml_setSpecies(SMOD, fba)
    sbml_setCompartments(SMOD, C)
    sbml_setReactions(SMOD, fba, return_dict=False)
    SMOD.setAnnotation(L3FBA.replace('<?xml version="1.0" encoding="UTF-8"?>',''))
    SDOC = sbml_writeSBMLModel(SMOD, fname, directory=directory, return_doc=True)

    if __DEBUG__:
        print '\nModel'        
        print L3FBA
        
        import pprint
        PP = pprint.PrettyPrinter()
        
        print '\nUnits'
        PP.pprint(sbml_setUnits(SMOD, units=None, give_default=True))
        print '\nCompartments'
        c,s = sbml_setSpecies(SMOD, fba, return_dicts=True)
        PP.pprint(c)
        print '\nSpecies'        
        PP.pprint(s)
        print '\nReactions'
        PP.pprint(sbml_setReactions(SMOD, fba, return_dict=True))

def writeSBML2FBA(fba, fname, directory=None, sbml_level_version=None):
    """
    Takes an FBA model object and writes it to file as SBML L3 FBA:

     - *fba* an fba model object
     - *fname* the model will be written as XML to *fname*
     - *directory* [default=None] if defined it is prepended to fname
     - *sbml_level_version* [default=None] a tuple containing the SBML level and version e.g. (2,4)

    """
    sbml_writeSBML2FBA(fba, fname, directory, sbml_level_version)
    